 subroutine UHFM(Hcore,U,nbasis,S,Enu)
        implicit none
        integer :: ii,jj,kk,ll,i,j,k,l
        integer :: mm,nn,kocc,kcy,nbasis
        integer :: nbeta,nalpha,maxcycle
        integer :: info,info1,info2
        real*8  :: rmsold,rms,Enu
        real*8  :: Eele,Eold,Ehf,Ethreshold,rmsthreshold
        real*8  :: Hcore(nbasis,nbasis),S(nbasis,nbasis)
        real*8  :: X(nbasis,nbasis),U(nbasis,nbasis,nbasis,nbasis)
        real*8  :: guessa(nbasis,nbasis),guessb(nbasis,nbasis),&
                   focka(nbasis,nbasis),fockb(nbasis,nbasis),&
                   fock0a(nbasis,nbasis),fock0b(nbasis,nbasis),&
                   Palpha(nbasis,nbasis),Pbeta(nbasis,nbasis),&
                   cofalpha(nbasis,nbasis),cofbeta(nbasis,nbasis),&
                   cof0alpha(nbasis,nbasis),cof0beta(nbasis,nbasis),&
                   orbalpha(nbasis),orbbeta(nbasis),Ptol(nbasis,nbasis),&
                   Ptolold(nbasis,nbasis), Palphaold(nbasis,nbasis),Pbetaold(nbasis,nbasis)
        real*8 ::tem(nbasis,nbasis),sem12(nbasis,nbasis),temmm1(nbasis,nbasis)
!!!!! for dsyev and dsyevd                
        real*8 ::work(3*nbasis+1),work1(1+6*nbasis+2*nbasis**2),work2(3+5*nbasis),&
                    eigv(nbasis)
!!!!
        maxcycle=10000
        nbeta=7
        nalpha=9
        Ethreshold=10**(-12.d0)
        rmsthreshold=10**(-12d0)
        focka=0d0
        fockb=0d0
        rmsold=0d0
        guessa=0d0
        guessb=0d0
!!!!  Diagonalize S,X=S^(-1/2)
!!!!!!!  dsyev is from lapack.
        Tem=S
        call dsyev('V','L',nbasis,Tem,nbasis,eigv,work,3*nbasis+1,info1)  
        ! call dsyevd('V','L',nbasis,Tem,nbasis,eigv,work1,1+6*nbasis+2*nbasis**2,work2,3+5*nbasis,info1)  
        if(info1/=0)stop " Diagonalize S failed !!!!! "
        Sem12=0d0
        do i=1,nbasis
            Sem12(i,i)=eigv(i)**(-0.5d0)
        end do   
!!!! Built-in functions  
        ! X=matmul(tem,Sem12)
        ! X=matmul(X,transpose(tem))
!!!! BLAS        
        call dgemm('N','N',nbasis,nbasis,nbasis,1d0,tem,nbasis,Sem12,nbasis,0d0,temmm1,nbasis)
        call dgemm('N','T',nbasis,nbasis,nbasis,1d0,temmm1,nbasis,tem,nbasis,0d0,X,nbasis)
!!!!  Build the Initial Guess Density
!!! for alpha:core Hamiltonian   
            guessa= Hcore   
!!! for beta:generalized Wolfsberg−Helmholz (GWH) 
!!! guess(i,j)=0.5*K*(Hcore(i,i)+Hcore(j,j))*S(i,j),K=1.75 typically
            do i=1,nbasis
                do j=1,i
                    if(i==j)then
                        guessb(i,j)=hcore(i,i)
                    else
                       guessb(i,j)=0.5d0*1.75*S(i,j)*(Hcore(i,i)+Hcore(j,j)) 
                       guessb(j,i)=guessb(i,j)
                    end if    
                end do
            end do
        focka=guessa
        fockb=guessb
        temmm1=matmul(transpose(X),Focka)
        fock0a=matmul(temmm1,X)
        temmm1=matmul(transpose(X),Fockb)
        fock0b=matmul(temmm1,X)
!!!!  
        tem=Fock0a
        call dsyevd('V','L',nbasis,Tem,nbasis,eigv,work1,1+6*nbasis+2*nbasis**2,work2,3+5*nbasis,info2)  
        if(info2/=0)stop " Diagonalize F alpha failed !!!!! "     
        cof0alpha=matmul(X,Tem)
        tem=Fock0b
        call dsyevd('V','L',nbasis,Tem,nbasis,eigv,work1,1+6*nbasis+2*nbasis**2,work2,3+5*nbasis,info2)  
        if(info2/=0)stop " Diagonalize F beta failed !!!!! "     
        cof0beta=matmul(X,Tem)
!!!!
!!! Density martrix  (alpha)
        cofalpha=transpose(cof0alpha)
        Palpha=0d0
          do mm=1,nbasis
              do nn=1,nbasis
                  do kocc=1,nalpha
                      Palpha(mm,nn)=Palpha(mm,nn)+cof0alpha(mm,kocc)*cofalpha(nn,kocc)
                  end do
              end do
          end do 
        !Palpha=matmul(cof0alpha(:,1:nalpha),transpose(cof0alpha(:,1:nalpha)))  
!!! Density martrix  (beta)
        cofbeta=transpose(cof0beta)
        Pbeta=0d0
        do mm=1,nbasis
            do nn=1,nbasis
                do kocc=1,nbeta
                    Pbeta(mm,nn)=Pbeta(mm,nn)+cof0beta(mm,kocc)*cofbeta(nn,kocc)
                end do
            end do
        end do         
        !Pbeta=matmul(cof0beta(:,1:nbeta),transpose(cof0beta(:,1:nbeta)))          
        Ptol=Palpha+Pbeta
!!!!        
        Eele=0d0
        do mm=1,nbasis
            do nn=1,nbasis
                Eele=Eele+Ptol(mm,nn)*Hcore(mm,nn)+Palpha(mm,nn)*focka(mm,nn)+Pbeta(mm,nn)*fockb(mm,nn)
            end do
        end do 
!!!! Start SCF    
        do kcy=1,maxcycle
            write(*,*)'iteration num :',kcy
            Eold=Eele
            Ptolold=Ptol
            Palphaold=Palpha
            Pbetaold=Pbeta
            rmsold=rms
            focka=0d0
            fockb=0d0
            do ii=1,nbasis
                do jj=1,nbasis
                    Focka(ii,jj)=Hcore(ii,jj)
                    Fockb(ii,jj)=Hcore(ii,jj)
                    do kk=1,nbasis
                        do ll=1,nbasis
                        ! ikjl=indexeri(ii,ll,kk,jj)
                        Focka(ii,jj)=Focka(ii,jj)+Ptol(kk,ll)*U(ii,jj,kk,ll)-Palpha(kk,ll)*U(ii,kk,jj,ll) 
                        Fockb(ii,jj)=Fockb(ii,jj)+Ptol(kk,ll)*U(ii,jj,kk,ll)-Pbeta(kk,ll)*U(ii,kk,jj,ll)                  
                        end do
                    end do
                end do
            end do   
!!!!!!! Diagonalize F'
            fock0a=matmul(matmul(transpose(X),focka),X)
            fock0b=matmul(matmul(transpose(X),fockb),X)
            tem=Fock0a
            orbalpha=0d0
            call dsyevd('V','L',nbasis,Tem,nbasis,orbalpha,work1,1+6*nbasis+2*nbasis**2,work2,3+5*nbasis,info)  
            if(info/=0)stop " Diagonalize F' failed !!!!! " 
            cofalpha=matmul(X,tem)
            tem=Fock0b
            orbbeta=0d0
            call dsyevd('V','L',nbasis,Tem,nbasis,orbbeta,work1,1+6*nbasis+2*nbasis**2,work2,3+5*nbasis,info)  
            if(info/=0)stop " Diagonalize F' failed !!!!! " 
            cofbeta=matmul(X,tem)            
!!!!!!! build new Density matrix           
!!! Density martrix  (alpha)
            Palpha=matmul(cofalpha(:,1:nalpha),transpose(cofalpha(:,1:nalpha)))          
!!! Density martrix  (beta)
            Pbeta=matmul(cofbeta(:,1:nbeta),transpose(cofbeta(:,1:nbeta)))   
            Ptol=Palpha+Pbeta
!!!! damping
            ! Ptol=(1d0-damp)*Ptol+damp*Ptolold  
            ! Palpha=(1d0-damp)*Palpha+damp*Palphaold  
            ! Pbeta=(1d0-damp)*Pbeta+damp*Pbetaold
!!!!            
            Eele=0d0
            rms=0d0
            do mm=1,nbasis
                do nn=1,nbasis
                    Eele=Eele+0.5d0*(Ptol(mm,nn)*Hcore(mm,nn)+Palpha(mm,nn)*focka(mm,nn)+Pbeta(mm,nn)*fockb(mm,nn))
                    rms=rms+sqrt((Ptol(mm,nn)-Ptolold(mm,nn))**2)
                end do
            end do 
            rms=rms/nbasis
                write(*,*)Eele-Eold, rms
            if(abs(Eele-Eold)<=Ethreshold .and.abs(rms)<=rmsthreshold)then
                write(*,*)'ok'
                exit
            end if
        end do !!! cycle
!!!!!!!!!!!!!  Outoput some information        
        Ehf=Enu+Eele
           write(*,*)' Final Hartree-Fock Energy :',Ehf,'A.U.'
    end subroutine UHFM
